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A  Monte  Carlo  Simulation  Approach 


FOREWORD 


A  technique  is  described  for  preserving  in  cost  estimates 
more  of  the  information  available  to  the  engineers  and  cost 
analysts  participating  in  the  process  of  estimating.  The  added 
information,  principally  in  the  form  of  a  frequency  distribution 
for  system  costs,  should  prove  useful  to  the  persons  faced 
with  the  problem  of  using  the  cost  estimate  in  the  decision¬ 
making  process. 

An  actual  test  of  the  technique,  made  with  a  major  com¬ 
ponent  of  the  Main  Battle  Tank  70,  will  be  documented  in  a  sub¬ 
sequent  publication. 


Arnold  Proschan 
Head,  Economics  and  Costing  Department 
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A  Monte  Carlo  Simulation  Approach 
to  Cost- Uncertainty  Analysis 


ABSTRACT 


An  Important  aspect  of  cost  research  is  the  measurement  of  the  un¬ 
certainty  Inherent  in  the  projection  of  system  cost.  Approaches  to  this 
problem  have  In  the  pastfocused  on  the  decisionmaker’s  intuition  oron  sen¬ 
sitivity  analysis.  Only  recently  have  approaches  utilizing  such  tools  as 
statistical  decision  theory  and  probability  theory  been  formulated. 

This  study  focuses  on  die  Monte  Carlo  simulation  approach  to  uncer¬ 
tainty  In  cost  analysis.  This  approach  requires  (a)  expression  of  input 
estimates  as  probability  distributions  reflecting  uncertainty  and  (b)  cost 
equations  pertinent  to  a  particular  model. 

The  Monte  Carlo  simulation  approach  then  generates  (a)  the  fre¬ 
quency  distribution  for  system  cost  and  (b)  statistical  measures  that  illus¬ 
trate  the  nature  and  magnitude  of  system  cost  uncertainty. 

Two  models  are  developed,  the  Beta  model  and  the  Weibull  model, 
each  of  which  reflects  a  particular  distribution  form  for  the  Inputs.  The 
relative  costs  and  advantages  of  each  model  are  compared. 

A  user’s  guide  to  the  program  and  complete  program  listings  are 
presented  In  App  A. 


INTRODUCTION 


The  purpose  of  this  paper  is  to  describe  a  technique  for  quantifying  un¬ 
certainty  in  cost  analysis.  Approaches  to  the  measurement  of  the  uncertainty 
present  in  cost  models  have  in  the  past  generally  centered  on  the  intuition  of 
the  decision  maker  or,  at  best,  on  sensitivity  analysis.  Only  recently  have  ap¬ 
proaches  based  on  statistical  decision  theory  and  probability  theory  been  at¬ 
tempted. 

The  technique  of  uncertainty  analysis  is  defined  for  purposes  of  this  paper 
as  the  use  of  probability  distributions  as  inputs  to  aggregate  models.  This 
analysis  can  be  used  in  several  ways.  For  example,  it  can  be  valuable  where 
there  is  a  major  lack  of  Information  concerning  the  value  of  input  parameters. 
This  can  be  true  of  systems  not  yet  in  existence  or  of  estimates  of  historical 
parameters  for  which  data  are  fragmentary.  It  can  also  be  of  value  in  esti¬ 
mating  the  cost  of  an  existing  system  for  which  actual  data  are  available. 
Samples  of  these  data  can  be  used  to  construct  distributions  for  input  variables. 

The  model  developed  is  called  the  “Monte  Carlo  Uncertainty  Analysis 
Model.”  It  is  to  be  used  with  any  cost  model.  The  model  allows  the  user  to 
specify  probability  distributions  for  the  cost  model  input  variables  rather  than 
the  usual  single-point  estimates.  These  are  then  input  to  the  cost  model  in  the 
uncertainty  analysis  model,  and  the  outputs  are  probability  distributions  and 
confidence  intervals  rather  than  single-point  estimates  of  costs.  There  are 
two  variations  of  the  model,  differing  as  to  types  of  probability  distributions 
that  can  be  used  to  characterize  inputs.  One  assumes  a  Beta  distribution  and 
the  other  a  Weibull,  The  availability  of  two  types  of  distribution  permits  the 
user  more  flexibility  in  specification.  The  rationale  for  use  of  one  or  the  other 
approach  rests  in  which  type  of  distribution  the  analyst  feels  most  closely  ap¬ 
proximates  the  distribution  of  the  input  variables  of  the  particular  problem  and 
which  set  of  input  specifications  he  can  most  easily  and  accurately  provide. 

The  first  chapter  of  the  paper  discusses  the  value  of  uncertainty  analysis 
models  in  general  and  describes  the  basic  methodology  of  techniques  developed. 
The  second  chapter  discusses  actual  use  of  the  techniques  in  several  simpli¬ 
fied  examples. 

The  appendix  provides  the  user  with  the  information  necessary  to  effec¬ 
tively  implement  the  programs.  The  precise  format  of  user-provided  inputs 
is  described.  Error  messages  are  incorporated  into  the  program  to  ensure 
conformity  with  input  requirements.  A  description  of  program  logic  is  fol¬ 
lowed  by  a  complete  listing  of  the  Weibull  and  Beta  programs. 


Chapter  1 


COST -UNCERTAINTY  ANALYSIS:  AN  OVERVIEW 


COST  UNCERTAINTY  IN  DECISION  MAKING 

The  process  of  decision  making  and  thus  the  fundamental  task  of  the  de¬ 
cision  maker  is  to  choose  among  alternative  courses  of  action.  Often  such 
choice  will  involve  a  cost- benefit  analysis  of  either  a  formal  or  an  informal 
nature.  Within  the  context  of  a  formal  analysis,  what  has  been  given  the  deci¬ 
sion  maker  as  costs  may  now  be  examined. 

The  precise  calculation  of  costs  is  not  a  difficult  task.  The  availability 
of  computers  and  cost  models  can  reduce  this  to  routine.  The  calculation  it¬ 
self  is  precise  and  can  be  done  rapidly  in  minute  detail. 

The  inputs  to  these  precise  calculations  are,  however,  not  precise;  in 
fact  they  are  often  quite  the  opposite.  The  errors  present  in  each  input  are 
passed  on  to’various  aggregations  until  a  total  cost  is  arrived  at  that  somehow 
reflects  each  individual  error.  Thus,  cost  data  inputs  are  combined  in  a  com¬ 
puter  model  with  a  multiplicity  of  equations  and  hundreds  of  other  inputs  to 
form  a  single  estimate  of  total  costs.  This  single  estimate  is  presented  to  the 
decision  maker  with  the  implication  that,  although  it  may  not  be  perfect,  it  is 
certainly  the  best  available  estimate  and  without  any  statement  as  to  likelihood 
of  occurrence  or  range  of  other  possible  values. 

But  in  the  generation  of  this  aggregate  number,  hundreds  of  imprecise 
numbers  may  have  been  used.  This  fact  is  not  emphasized  to  the  decision 
maker  nor  does  he  have  a  basis  on  which  to  judge  the  precision  of  the  numbers. 
Cost  models  as  currently  conceived  and  used,  then,  uniformly  withhold  from 
the  decision  maker  some  information  that  might  be  vital  to  his  decision.  He 
has  been  denied  some  available  information  on  the  accuracy  of  the  estimates. 

As  an  example,  suppose  that  one  of  the  elements  to  be  costed  is  a  missile 
airframe.  The  missile  is  not  yet  designed  or  built,  so  no  production-cost  data 
are  available.  A  cost-estimating  relation  (CER)  could  be  used,  for  example, 
to  estimate  that  the  airframe  cost  for  a  new  missile  would  be  about  $32,000. 
The  design  group,  however,  warns  that  some  added  sophistication  might  make 
the  CER  predict  on  the  low  side,  but  the  improvements  in  some  manufacturing 
techniques  promise  lower  costs.  Some  quick  calculations  show  estimates  as 
low  as  $26,000  and  a3  high  as  $45,000,  Each  of  these  calculations  is  based  on 
a  set  of  assumptions  concerning  labor  and  overhead  rates,  material  costs  in 
the  future,  and  design  details  not  yet  firm  and  subject  to  some  uncertainty. 


For  each  of  a  multitude  of  other  important  inputs,  similar  assumptions 
are  made,  and  a  single  aggregate  cost  is  obtained  as  the  output  of  the  cost 
model.  The  hypothetical  decision  maker  is  then  forced  to  make  his  final  de¬ 
cision  on  a  single  value  derived  from  several  imprecise  input  values.  If  he 
can  be  provided  with  some  additional  knowledge  concerning  the  impact  of  the 
uncertain  variable,  a  more  rational  decision  can  be  made. 

The  dangers  of  the  single-cost  estimate  have  been  recognized  for  years, 
and  several  strategies  have  been  developed  for  augmenting  the  analysis  or 
circumventing  the  difficulty.  One  is  isolation  of  the  differences  between  alter¬ 
natives.  Cost  elements  common  to  alternatives  are  estimated  in  a  similar  or 
identical  manner  so  that  only  the  uncertainties  of  the  unique  feature  of  alter¬ 
natives  affect  relative  cost.  Another  is  the  use  of  sensitivity  analysis.  In 
sensitivity  analysis  the  impact  of  errors  in  estimates  and  assumptions  is  com¬ 
puted,  and  error  sources  important  to  the  choice  of  alternatives  are  identified. 
Such  an  analysis  can  produce  proof  in  insensitivity,  or  it  can  provide  evidence 
that,  within  a  relevant  range  of  values,  choice  is  or  is  not  affected  by  estima¬ 
tion  error. 

When  judgments  about  the  relevant  range  of  all  variables  in  a  model  are 
considered,  sensitivity  analysis  produces  an  array  of  numbers  that  includes 
the  analysts’  beliefs  concerning  the  limits  of  the  variables  but  excludes  any 
knowledge  of  their  relative  probability.  No  probability  statements  are  fur¬ 
nished,  and  the  decision  maker  could  be  led  to  believe  that  all  numbers  in  the 
array  are  equally  likely. 

The  choice  of  a  relevant  range  of  values  for  sensitivity  analysis  is  both 
difficult  and  critical  to  its  usefulness.  It  also  reveals  a  dilemma  inherent  in 
the  application  of  sensitivity  analysis.  If  the  analyst  has  evidence  that  the 
value  of  one  of  the  inputs  is  constrained  within  some  upper  and  lower  limits, 
then  this  same  evidence  may  provide  information  on  the  relative  likelihood  of 
particular  values.  In  some  cases  an  intuitive  belief  about  the  range  may  pro¬ 
vide  him  with  equally  valid  suppositions  concerning  probabilities  of  the  values 
within  the  range. 

In  the  previous  example  of  the  missile  cost,  if  the  reasoning  that  fixed 
the  relevant  range  of  costs  at  $26,000  to  $45,000  could  yield  information  on 
the  probability  of  occurrence,  the  missing  information  could  be  supplied  in  a 
form  useful  to  the  decision.  A  possible  statement  might  be  that  there  is  a  99 
percent  probability  that  the  cost  is  less  than  $42,000. 

When  all  uncertain  elements  of  a  given  cost  model  are  combined,  sensi¬ 
tivity  analysis  may,  for  example,  produce  a  maximum  cost  of  $810  million 
even  though  there  is  only  one  chance  in  a  hundred  that  costs  will  exceed  $700 
million.  In  the  sensitivity  analyses,  if  the  $800  million  number  exceeded 
slightly  the  cost  of  another  alternative  under  similar  sensitivity  assumptions, 
the  decision  maker  may  have  been  furnished  an  unlikely  cost  set  for  considera¬ 
tion  along  with  all  other  sensitivity  sets,  and  this  set  may  uniquely  favor  a 
different  alternative.  Knowing  just  a  few  simple  probability  statements,  then, 
can  make  unnecessary  the  time-consuming  consideration  of  cost  estimates 
whose  likelihood  is  very  remote. 

In  short,  although  sensitivity  analysis  is  a  powerful  tool  for  portraying 
the  results  of  estimating  error,  it  may  lead  to  considering  highly  unlikely 


situations.  The  same  reasoning  that  leads  to  a  determination  of  relevant 
range  for  sensitivity  analysis  has  the  potential  of  providing  key  information 
in  decision  making. 

Value  of  Probability  Information 

The  use  of  probability  distributions  as  inputs  to  aggregate  models  has 
recently  captured  the  attention  of  the  cost  analysts.  The  two  most  widely  used 
techniques  for  handling  subjective  probability  information  are  the  derivation 
of  moments1  and  the  Monte  Carlo  simulation.  This  paper  focuses  on  the  Monte 
Carlo  simulation  approach  and  relies  heavily  on  the  works  of  David  B.  Hertz,2 
Paul  F.  Dienemann,3  and  W.  D.  Lamb.4 

A  discussion  of  the  Monte  Carlo  m  lology  is  preceded  by  several  ex¬ 
amples  designed  to  illustrate  the  potential  value  of  probability  distribution  in¬ 
formation  to  the  decision  maker.  In  each  illustration  the  frequency  distribu¬ 
tions  for  two  alternatives  are  shown.  The  horizontal  axis  in  each  case  repre¬ 
sents  the  cost  of  the  alternative  and  is  increasing  to  the  right.  The  vertical 
axis  represents  the  likelihood  of  occurrence  at  each  cost  level.  Each  of  these 
is  a  hypothetical  case  in  which  equal  effectiveness  or  other  benefits  are  as¬ 
sumed.  The  decision  maker’s  problem  is  that  of  choosing  the  least-cost 
alternative.  If  only  single-point  cost  estimates  were  provided  the  decision 
maker  would  of  course  feel  constrained  to  select  the  lower  cost  in  each  case. 
However,  Fig.  1  demonstrates  how  information  provided  by  probability  estima¬ 
tion  could  modify  his  outlook. 

The  peak  of  each  curve  is  at  the  most  likely,  or  modal,  value  and  that  is, 
in  these  hypothetical  cases,  the  only  cost  total  that  would  be  furnished  a  deci¬ 
sion  maker  in  the  absence  of  uncertainty  analysis.  In  thq  example  in  Fig.  1 
B  is  expected  to  be  less  expensive  but  it  has  a  much  larger  variance,  so  that 
extremely  high  costs  are  more  likely  than  with  A.  Faced  with  this  dilemma, 
the  decision  maker  may  decide  to  avoid  extreme  costs  by  choosing  A  or  to 
gamble  on  the  expected  lower  costs  of  B.  The  decision  cannot  be  prescribed 
for  him,  but  probability  distributions  can  give  him  awareness  of  a  pitfall  in 
selecting  the  “less  expensive”  alternative. 

Figure  2  illustrates  a  clear-cut  case.  The  largest  possible  cost  of  one 
alternative  is  less  than  the  smallest  possible  cost  of  the  other.  Decision 
makers  furnished  with  this  information  are  not  likely  to  choose  differently 
than  if  only  point  estimates  were  given. 

In  Fig.  3  a  single-point  estimate  would  furnish  no  basis  for  choice.  With 
most  likely  values  essentially  the  same  the  decision  maker  must  look  else¬ 
where  for  t  ifferences.  If  the  probability  distributions  are  furnished,  however, 
it  is  apparent  that  the  costs  cannot  be  regarded  as  equal.  If  they  appear  as 
shown  here  it  is  presumed  that  A  would  be  chosen,  since  it  can  be  better  ac¬ 
commodated  in  the  budgeting  and  financial  management  system.  It  is  possible 
that  much-lower -than-mean  costs  may  have  some  value  too,  and  as  a  result 
it  becomes  impossible  to  speculate  on  choice  outside  specific  cases.  The  im¬ 
portant  thing  is  that  information  helpful  in  the  decision  process  has  been  fur¬ 
nished  that  would  not  otherwise  be  available. 

The  costs  associated  with  two  alternatives  may  or  may  not  be  signifi¬ 
cantly  different.  Figure  4  illustrates  two  cost  totals,  each  with  the  same 
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variance.  If  variance  is  low  A'  and  B',  the  situation  is  as  in  Case  2;  costs  are 
different  and  offer  a  basis  for  choice  as  conclusive  as  the  magnitude  of  the 
difference  might  indicate.  Larger  variances,  A  and  B,  diminish  the  importance 
of  the  most-likely -value  cost  difference.  At  some  combination  of  “closeness* 
and  high  variance  the  cost  difference  may  not  be  of  significance  in  the  selection 
of  an  alternative.  In  any  case,  quantification  of  the  probabilities  of  the  differ¬ 
ences  is  a  useful  contribution  to  the  decision  maker’s  understanding. 

These  four  examples,  it  should  be  noted,  differ  only  quantitatively.  They 
are  offered  in  these  forms  to  illustrate  the  range  of  possibilities  in  which  useful 
information  may  come  from  a  knowledge  of  the  probability  distributions  of 
total  cost. 


I  METHODOLOGY 

!? 

i 

!  This  section  presents  the  basic  methodology  of  two  variations  of  the  Monte 

|  Carlo  model.  One  assumes  that  uncertain  inputs  can  be  characterized  by  Beta 

{  distributions.  The  other  variation  uses  the  Weibull  distribution. 


Overview  of  Logic 

Since  the  two  variations  of  the  Monte  Carlo  technique  are  identical  as  to 
logic  and  statistical  methods,  the  general  logic  of  both  models  can  be  discussed 
together.  Basically  the  models  consist  of  a  Monte  Carlo  simulation,  a  listing 
of  the  simulation  results  by  class  intervals,  and  a  plot  of  the  output  distribu¬ 
tions.  In  the  simulation  the  value  of  each  input  parameter  is.  chosen  randomly, 
based  on  the  probability  distribution  of  the  parameter.  When  this  procedure 
is  repeated  many  times,  a  range  of  values  is  produced.  These  values  are  aggre¬ 
gated  into  a  number  of  class  intervals  and  results  are  plotted.  Finally,  since 
the  results  of  the  simulation  produce  a  very  crude  plot,  a  smoothing  routine 
produces  a  more  regular  version  of  the  output  distribution. 

Beta  Variant 

This  approach  is  quite  similar  to  one  formulated  by  P.  Dienemann  of  The 
RAND  Corporation.3  It  assumes  that  cost  distributions  can  be  approximated 
by  Beta  distributions.  These  are  unimodal  and  continuous  finite  at  both  upper 
and  lower  bounds  and  can  be  skewed  or  symmetric.*’5  This  section  discusses 
(a)  required  inputs,  (b)  general  logic  of  each  technique,  and  (c)  outputs. 

Inputs.  To  specify  the  Beta  distribution  for.  an  input  variable  the  following 
input  parameters  are  required: 

(1)  XP  (the  most  probable  value  of  the  distribution!:  If  one  views  the 
probability  distribution  of  an  X,  Y  plane,  the  most  probaMo  value  corresponds 
to  the  point  having  the  greatest  Y  value.  For  example,  the  points  P  and  P,  are 
the  most  probable  values  for  the  distributions  shown  in  Fig.  5.  In  the  figure 
the  abscissa  describes  the  values  that  the  input  variable  may  assume  and  the 


*A  more  comprehensive  discussion  of  this  type  of  distribution  can  be  found  in 
Mood  and  Grayblll,®  p  129. 
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ordinate  gives  the  likelihood  of  occurrence  of  each  of  those  values.  In  the 
symmetric  distribution  (curve  A  of  Fig.  5)  the  most  probable,  or  modal,  value 
as  it  is  generally  known  also  corresponds  to  the  expected  or  mean  value.* 

(2)  XH  (high  value):  This  is  the  point  to  the  right  of  the  modal  value 
where  the  probability  distribution  meets  the  abscissa  given  the  X  values  in¬ 
crease  from  left  to  right.  In  Fig.  5  the  points  H  and  Hi  are  the  high  values. 


Fig.  5 — Typical  Beta  Distributions 
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Fig.  6 — Sample  Set  of  Distributions 

(3)  XL  (low  value):  This  is  the  point  to  the  left  of  the  modal  value  where 
the  probability  distribution  meets  the  abscissa  given  the  X  values  increase 
from  left  to  right.  In  Fig.  5  the  points  L  and  L:  are  the  low  values.  In  Fig.  5 
they  are  the  zero  point  of  a  standard  X,  Y  coordinate. 

(4)  Distribution  Type:  Finally  the  user  must  compare  his  concept  of  the 
probability  distribution  of  the  input  variable  with  a  set  of  standard  distributions 
and  choose  the  one  that  most  closely  approximates  his  concept.  The  set  o' 
distributions  now  being  used  is  shown  in  Fig.  6,  These  distributions  can  be 
described  qualitatively  on  the  basis  of  (a)  being  symmetric  or  skewed  or  (b) 
having  different  degrees  of  variances.  This  is  illustrated  in  Fig.  6.  Most  un- 

*For  a  strict  definition  of  this  see  Mood  and  Graybill,5  p  103. 


certain  inputs  should  be  approximated  reasonably  by  one  of  the  nine  distribu¬ 
tion  types. 

(5)  The  Cost  Model:  The  user  must  provide  the  costing  equations  that 
compose  his  cost  model.  These  are  embedded  within  the  uncertainty  analysis 
model  and  are  necessary  to  compute  the  values  from  which  the  final  plot  of  the 
distributions  is  to  be  made.  To  illustrate  with  an  extremely  simple  example, 
suppose  one  were  interested  in  calculating  the  probability  distribution  of  per¬ 
sonnel  and  maintenance  costs  for  an  aircraft  battalion.  Further,  suppose  per¬ 
sonnel  costs  were  available  on  a  per  capita  basis  and  maintenance  on  a  per 
aircraft  basis,  then  the  cost  model  submitted  by  the  user  might  look  like  this: 

Fl  =  number  of  men  x  dollars  per  man 

F2  =  number  of  aircraft  x  maintenance  cost  per  aircraft 

Total  cost  =  Fl  +  F2 

The  five  types  of  input  described  above  must  be  prepared  for  each  input 
variable  In  the  cost  model.  The  necessary  model  inputs  having  been  discussed, 
the  general  logic  of  the  model  may  now  be  considered. 


For  various 
levels  of 


Fig.  7— Monte  Carlo  Technique,  Beta  Variant 


Logic  and  Output.  Figure  7  illustrates  the  basic  methodology  of  the  model. 
For  each  of  the  uncertain  costs  inputs,  the  four  required  parameters  are  pro¬ 
vided  by  the  user.  Using  this  information,  the  model  randomly  selects  a  single 
value  for  each  input.  This  particular  set  of  values  is  then  input  to  the  cost  model 
that  is  now  embedded  in  the  uncertainty  model.  This  produces  a  particular  set 
of  output  values  that  are  stored  for  use  later  in  building  the  final  distribution 
curves.  The  process  is  repeated  1000  times,  which  should  be  a  sufficient 
number  of  iterations  to  develop  an  accurate  total  cost  distribution.  The  number 
of  iterations  can  be  changed  easily  to  meet  user  requirements.  If  the  specified 
number  of  iterations  has  not  been  run,  the  model  will  recycle  to  again  (a)  ran¬ 
domly  select  different  values  for  each  input,  (b)  compute  output  values  for  these 
selected  inputs  using  the  embedded  cost  equations,  and  (c)  store  these  results 
for  further  use.  ' 

When  the  specified  number  of  Iterations  has  been  run,  the  model  moves 
to  a  new  phase.  This  Involves  compilation  of  the  results  of  e  ch  run  to  (a) 
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prepare  frequency  distributions  for  aggregate  cost  and  (b)  compute  statistical 
measures  such  as  mean  and  variance.  First,  the  values  are  ordered  into  class 
intervals  covering  the  range  of  possible  values.  The  frequency  diagram  pro¬ 
duced  by  this  ordering  is  then  smoothed  to  produce  the  more  familiar  form  of 
probability  distribution  illustrated  in  Fig.  7.  The  next  chapter  contains  a  more 
detailed  description  of  the  actual  form  of  output  generated  by  the  model. 

Weibull  Variant 

In  this  approach  inputs  are  characterized  by  a  Weibull  distribution.*  >®>T 
Weibull  distributions,  like  Beta  distributions,  are  unimodal  continuous  and  can 
be  symmetric  or  skewed  but,  unlike  Beta  distributions,  are  infinite  at  both 
upper  and  lower  bounds.  Since  the  Weibull  variant  approach  is  identical  in 
logic  and  output  to  the  Beta  variant  just  described,  only  the  input  requirements 
of  the  Weibull  variant  will  be  discussed  in  detail. 


For  various 
levels  of 


Fig.  8 — Monte  Carlo  Technique,  Weibull  Variant 

Inputs.  Figure  8  illustrates  the  methodology  of  the  Weibull  variant.  As 
illustrated  in  the  figure  the  input  requirements  for  each  cost  input  are: 

(1)  XP  (most  probable  value):  This  has  the  same  meaning  as  the  most 
probable  value  input  for  the  Beta  variant.  The  discussion  is  identical  and 
therefore  will  not  be  repeated  at  this  point. 

(2)  XH  and  PH:  Since  these  inputs  are  interrelated,  both  can  be  dis¬ 
cussed  in  one  section.  XH  is  a  value  of  the  input  that  is  higher,  i.e.,  greater 
than,  the  most  probable  value  specified  for  the  input.  PIT  is  the  probability 
that  the  value  of  that  input  will  be  some  quantity  greater  than  XH.  Figure  9 
illustrates  this  graphically.  The  X  axis  is  the  value  of  the  input  variable  and 
the  Y  axis  is  the  probability  of  occurrence  of  different  X  values.  Upon  closer 
inspection,  it  is  clear  that  this  input  specification  is  analogous  to  the  specifica¬ 
tion  of  XH  for  the  Beta  variant.  In  the  case  of  the  Beta  variant  the  user  was 
asked  to  specify  XH  such  that  the  probability  of  the  input  being  greater  than 

*A  complete  mathematical  description  of  the  Weibull  distribution  can  be  found  in 
Refs  6  and  7. 
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XH  was  zero.  Thus  PH  is  zero.  Since  the  Weibull  distribution  is  infinite  and 
thus  has  no  fixed  bound  the  user  is  asked  to  specify  XH  and  PH  in  such  a  way 
that  the  probability  of  the  value  being  greater  than  XH  is  PH. 

(3)  XL  and  PL:  These  are  analogous  to  XH  and  PH  where  XL  is  some 
value  less  than  XP  such  that  the  probability  of  the  input  variable  being  less 
than  XL  is  PL. 

Logic  and  Outputs.  With  these  inputs  the  Weibull  variant  operates  in 
essentially  the  same  fashion  as  the  Beta  variant.  Values  are  randomly  se¬ 
lected  for  each  of  the  cost  inputs,  and  these  are  processed  through  the  cost 
model  with  results  being  stored  for  development  of  the  final  distributions. 

After  the  specified  number  of  iterations  is  completed  the  model  orders  the 
results  of  each  iteration  and  develops  frequency  distributions,  which  are 
smoothed  and  then  finally  graphed. 

Model  Limitations.  The  techniques  developed  in  this  paper  provide  many 
advantages  over  more  conventional  forms  of  cost  analysis.  Uncertainty  about 
input  variables  can  be  explicitly  described.  The  simultaneous  interaction  of 
uncertainties  in  many  variables  can  be  assessed  and  then  graphically  displayed 


to  those  confronted  by  decision.  The  sensitivity  of  final  results  to  each  or  a 
set  of  inputs  can  be  efficiently  tested.  However,  it  would  be  at  best  naive  to 
claim  there  are  no  disadvantages  to  use  of  these  techniques.  Some  major  con¬ 
siderations  are: 

(a)  The  assumption  of  Beta  or  Weibull  distribution  may  not  be  a  sound  one. 

(b)  Much  additional  input  information  is  necessary  to  operate  the  model. 

(c)  The  model  assumes  independence  of  the  input  variables.  This  can 
prove  limiting  in  certain  situations.  There  are  at  present  several  approaches 
to  handling  dependency.  Among  these  are  incorporation  in  the  model  of  the 
functional  relation  between  the  variables,  statement  of  the  dependent  variable 
in  terms  of  auxiliary  variables,  and  use  of  joint  probability  distributions. 

Though  these  can  eliminate  dependence  in  a  majority  of  costing  problems,  the 
systematic  handling  of  input  dependency  remains  a  very  open  area  for  research. 


Chapter  2 

IMPLEMENTATION  OF  MONTE  CARLO  MODELS 


A  thorough  understanding  of  the  input  requirements  and  output  interpreta¬ 
tion  of  uncertainty  analysis  models  is  necessary  before  attempting  to  apply  the 
models  to  specific  problems.  The  user  of  the  models  should  be  aware  of  model 
limitations,  sources  of  input,  output  interpretation,  and  the  differences  between 
the  Beta  and  Weibull  models  in  order  to  implement  the  models  effectively.  A 
discussion  of  computer  time  requirements  for  the  models  is  presented  to  ac¬ 
quaint  the  potential  user  with  the  anticipated  costs  of  applying  the  program  to 
a  given  problem.  Format  of  inputs,  computer  program  logic,  and  complete 
program  listings  are  presented  in  the  appendix. 


SOURCES  OF  INPUT  INFORMATION 

The  Beta  and  Weibull  simulation  models  are  designed  to  handle  any  type 
of  cost  function.  Information  concerning  uncertain  variables  may  come  from  a 
variety  of  sources.  The  user  of  the  model  must  be  prepared  to  obtain  the  pre¬ 
cise  input  information  required  by  the  model  from  these  diverse  sources. 

In  many  cases  the  best  source  of  information  about  an  uncertain  input  is 
a  person  or  agency  familiar  with  the  variable  in  question  and  possible  sources 
of  variation.  For  example,  a  certain  cost  model  may  require  the  cost  of  a 
newly  designed  equipment  item  as  an  input.  An  engineer  who  is  working  on 
the  design  of  the  equipment  and  who  is  acquainted  with  the  possibilities  of  pro¬ 
duction  delays  and  future  design  modifications  may  be  the  best  source  for  nec¬ 
essary  inputs.  The  input  is  still  an  estimate,  but  by  minds  more  intimately 
informed  on  the  details  of  a  specific  variable. 

Another  example  is  a  cost  model  that  requires  manpower  level  as  an  in¬ 
put.  The  best  source  of  information  may  be  an  experienced  military  advisor 
who  is  familiar  with  manpower  requirements  under  varying  circumstances. 

In  both  cases,  the  expert  providing  the  information  may  be  unfamiliar  with 
statistics  and  the  precise  meaning  of  model  input  requirements  such  as  finite 
upper  and  lower  limits  or  assignment  of  probabilities. 

It  is  necessary  for  the  user  to  precisely  define  the  variable  that  must  be 
estimated  and  to  carefully  delineate  the  information  needed  from  such  expert 
sources  of  information.  An  excellent  example  of  techniques  for  obtaining  pro¬ 
gram  inputs  from  such  sources  is  presented  in  “A  Technique  for  Probability 
Assignment  in  Decision  Analysis”  by  W.  D.  Lamb.4 
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Information  concerning  uncertain  inputs  is  frequently  derived  from  sta¬ 
tistical  techniques  such  as  regression  and  correlation  analysis.  Examples  are 
(a)  cost-estimating  relations  that  express  the  cost  of  equipment  or  material  as 
a  function  of  design  or  performance  characteristics,  (b)  overhaul  rates  for 
military  vehicles,  and  (c)  the  ratio  of  draftees  to  total  accessions  to  the  armed 
services.  Information  obtained  from  regression  models  can  be  readily  trans¬ 
lated  into  the  form  required  by  the  uncertainty  analysis  models.  Whether  the 
best  information  source  is  an  engineering  expert,  a  regression  model,  or  any 
combination  of  sources,  the  information  can  be  expressed  in  the  form  required 
for  an  input  in  the  uncertainty  model. 


MODEL  INPUT  REQUIREMENTS 
The  Beta  Model 

The  necessary  inputs  for  the  Beta  model  are: 

(1)  XP:  The  most  likely  value  of  the  variable 

(2)  XH:  A  finite  upper  limit 

(3)  XL:  A  finite  lower  limit 

(4)  Distribution  type* 

The  first  three  inputs  must  be  obtained  from  the  best  available  source; 
the  fourth  can  be  derived  from  the  first  three  by  the  user.  The  most  readily 
available  input  should  be  XP,  the  most  likely  value  of  the  variable.  Obtaining 
finite  upper  and  lower  limits  presents  a  more  serious  problem. 

In  many  cases  the  expert  source  of  information  will  hesitate  to  provide 
a  finite  upper  limit  for  an  uncertain  variable,  for  to  do  so  would  be  to  imply 
that  the  value  of  the  variable  cannot  possibly  exceed  this  limit.  Consider  the 
case  of  the  engineer  estimating  the  cost  of  an  advanced-design  engine  at  an 
early  stage  in  its  development.  The  remote  possibility  of  a  major  change  in 
end-item  design  characteristics  or  an  insoluble  technical  problem  could  drive 
the  cost  up  considerably.  These  possibilities  may  preclude  determination  of 
a  finite  upper  limit  or  require  an  extremely  high  upper  limit  to  account  for  any 
feasible  situation.  Similar  problems  could  be  encountered  by  the  military  ana¬ 
lyst  who  is  predicting  troop  strength  in  some  future  time  period.  The  prob¬ 
ability,  though  very  small,  of  an  armed  conflict  requiring  massive  troop  build¬ 
up  or  of  a  technological  breakthrough  that  will  replace  existing  manpower  re¬ 
quirements  dramatically  must  be  considered  in  establishing  the  limits  of  the 
distributloh. 

The  finite  nature  of  the  Beta  distribution  implies  that  the  value  of  the  un¬ 
certain  variable  must  remain  within  the  range  described  by  the  upper  and 
lower  limits  with  100  percent  certainty.  The  distribution  of  an  uncertain  var¬ 
iable  obtained  through  regression  analysis  with  normal  error  terms  has  in¬ 
finite  range.  The  model  user  can,  however,  use  a  99  percent  confidence  in¬ 
terval  as  a  reasonable  proxy  for  certainty. 

Once  the  most  likely,  high,  and  low  values  are  determined,  the  model  user 
can  derive  the  distribution  type.  If  the  difference  between  the  upper  limit  and 

^Distribution  types  characterized  by  direction  of  skewness  and  relative  variance 
are  presented  la  Fig.  8. 
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the  most  likely  value  (XH  -  XP)  equals  the  difference  between  the  most  likely 
value  anu  the  lower  limit  (XP  -  XL),  the  distribution  is  symmetric.  If  (XP  -  XL) 
is  greater  than  (XH  -  XP),  the  distribution  is  skewed  to  the  left.  If  (XH  -  XP) 
is  greater  than  (XP  -  XL),  the  distribution  is  skewed  to  the  right.  Once  the 
direction  of  skewness  is  determined,  the  user  can  ask  the  expert  source  to 
choose  from  the  three  figures  representing  the  relative  variance  of  distribu¬ 
tions  in  the  category.  If  the  information  is  obtained  from  a  regression  model, 
the  symmetric  medium  variance  curve  is  the  most  appropriate  since  it  approx¬ 
imates  the  transformation  of  a  normal  curve  to  a  Beta  distribution. 

The  Weibull  Model 

Weibull  input  requirements  are: 

(1)  XP:  The  most  likely  value 

(2)  XH:  The  upper  limit 

(3)  PH:  Probability  that  the  value  of  the  input  will  exceed  XH 

(4)  XL:  The  lower  limit 

(5)  PL:  Probability  that  the  value  of  the  input  will  be  less  than  XL 

Values  for  XP,  XL,  and  XH  are  obtained  as  in  the  preceding  model.  The 

major  difference  is  that  XL  and  XH  are  not  finite  limits.  Some  of  the  problems 
associated  with  the  Beta  model  are  thus  eliminated.  The  expert  source  may  be 
more  willing  and  able  to  present  bounds  with  a  probability  value  attached  to 
them.  When  uncertain  inputs  are  obtained  by  means  of  regression  analysis 
with  error  terms  normally  distributed,  the  Weibull  input  requirements  are 
easily  and  accurately  met.*  Since  the  variable  has  a  known  normal  distribu¬ 
tion,  XP  is  the  predicted  value  of  the  variable  and  PH  and  PL  are  chosen  by 
the  user  to  determine  XL  and  XH.  If  the  user  chose  5  percent  as  values  for 
both  PH  and  PL,  a  90  percent  confidence  interval  about  the  predicted  value 
would  yield  XL  and  XH. 


OUTPUT  FORMAT 

The  outputs  of  the  Beta  and  Weibull  Monte  Carlo  programs  are  identical 
in  format.  Both  routines  produce  two  plots  of  the  total  cost  distribution  and  a 
table  of  frequencies  within  class  intervals. 

The  first  plot  produced  by  the  program  shown  in  Fig.  10  is  that  of  the 
total  cost  distribution  as  produced  by  the  Weibull  Monte  Carlo  simulation. 

The  Y  axis,  running  horizontally  across  the  top  of  the  computer  output  page, 
measures  the  number  of  total  cost  estimates  lying  within  a  given  class  inter¬ 
val.  The  X  axis,  running  vertically  along  the  left-hand  side  of  the  computer 
output  page,  measures  total  cost.  The  “Xs”  that  form  the  distribution  show 
how  many  times  total  cost  fell  within  each  interval.  Four  other  values,  XMIN, 
XMEAN,  XMAX,  and  STD  DEV,  are  printed.  These  represent  the  minimum 
total  cost,  the  mean  total  cost,  the  maximum  total  cost,  and  the  standard  devia¬ 
tion  of  the  total  cost  distribution.  In  Fig.  10,  for  example,'  th'  minimum  cost 
is  $98,422,  the  mean  cost  is  $116,111,  the  maximum  cost  n3  $131,859,  and  the 
standard  deviation  is  $5697. 

*‘The  Weibull  density  function  is  nearly  symmetric  (approximating  the  normal  dis¬ 
tribution)  when  its  shape  parameter  is  approximately  3.5.”  See  Lamb,4  p  15. 
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X  AXIS  (COST  INTERVALS) 


The  second  plot,  Fig.  11,  is  a  smoothed  version  of  the  first.  The  cost 
data  generated  by  the  Monte  Carlo  simulation  are  smoothed  by  fitting  them  to 
a  polynomial  by  the  method  of  “least  squares.**  The  frequencies  are  rounded 
to  the  nearest  integer. 

The  table  of  frequencies  produced  by  the  program  and  shown  in  Fig.  12 
has  four  columns  of  information.  The  leftmost  column  contains  the  midpoint 
of  each  of  the  100  class  intervals.  Column  2  contains  the  number  of  total  cost 
estimates,  produced  by  the  Monte  Carlo  simulation  routine,  that  lie  within 
each  cost  interval.  The  third  column  contains  the  Y  value  produced  by  the 
smoothing  routine  for  each  class  interval.  This  value  is  rounded  to  the  near¬ 
est  integer.  The  fourth  column  is  the  difference  between  the  second  and  third 
columns.  In  Fig.  12,  for  example,  the  highest  class  interval  has  a  midpoint  of 
$131,692,  a  frequency  of  2.0,  a  smoothed  frequency  of  2.4291,  and  the  differ¬ 
ence  between  the  two  frequencies  is  0.4291. 


.  V 


COMPUTER  TIME  REQUIREMENTS 

The  purpose  of  this  section  is  to  provide  the  user  with  computer  time  re¬ 
quirements  as  a  measure  of  cost  to  enable  him  to  evaluate  these  programs  in 
terms  of  cost  effectiveness. 

For  the  Beta  and  the  Weibull  programs  two  parameters  were  varied  in 
order  to  determine  the  sensitivity  of  computer  time  consumption.  These 
parameters  are  the  number  of  iterations  in  the  Monte  Carlo  simulation  and 
the  number  of  variables.  It  was  considered  a  priori  that  these  two  parameters 
would  explain  most  of  the  time  variations  from  run  to  run. 

The  number  of  iterations  was  set  at  500  and  1000  for  the  purposes  of  ex¬ 
perimentation.  The  number  of  cost  variables  was  changed  from  3  to  15.  Along 
with  the  change  in  the  number  of  variables  there  was  a  change  in  the  complex¬ 
ity'  of  the  algebraic  operations.  These  two  cost  models  are  shown  in  Fig.  13. 

Timing  results  for  the  Beta  and  Weibull  programs  are  shown  in  Table  1. 
Total  time  results  include  execution  time,  compiler  time,  system  time,  and 
load  time.  Since  execution  time  varies  with  the  number  of  iterations  and  the 
compiler,  load,  and  system  times  are  roughly  fixed,  it  is  also  useful  to  see 
the  execution  times  alone  since  they  are  roughly  analogous  to  recurring  costs. 

A  number  of  points  should  be  noted  in  interpreting  Table  1.  First,  the 
times  shown  are  for  an  IBM  7044.  If  the  programs  are  used  on  a  different 
computer  the  times  will  vary.  Second,  the  total  times  are  quite  small,  never 
exceeding  5  min.  Third,  for  the  most  part  execution  time  is  less  than  50  per¬ 
cent  of  total  time.  Finally,  the  execution  time  will  not  vary  greatly  with  the 
number  of  variables  since  the  major  difference  is  the  number  of  data  cards 
that  must  be  read  by  the  computer. 

From  this  information  it  seems  probable  that  the  vast  majority  of  cost 
models  will  require  less  than  10  min  of  computer  time. 


*The  method  of  least  squares  chooses  the  coefficient  of  the  fitted  polynomial  so 
that  the  sum  of  the  squared  deviations  between  the  polynomial  and  real  data  Is  mini¬ 
mized. 
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Fig,  12 — Table  of  Frequencies  Plotted 
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Fig.  13— Cost  Models  a  and  b  Used  in  Generating  Timing  Requirements 


TABLE  1 


Computer  Time  Requirements 


Total  time,0  min 

Execution  time 

V  =3b 

V  =15 

V  =3 

1! 

> 

Beta 

Ic  =  500 

2.61 

3.07 

0.64 

0.87 

[  =1000 

3.17 

3.30 

0.74 

1.21 

ffeibull 

I  =500 

2.39 

2.67 

0.25 

0.44 

I  =1000 

2.49 

2.92 

0.35 

0.69 

aTotal  time  includes  compiler  time,  system  time,  load  time, 
and  execution  time.  All  times  are  measured  on  an  IBM  7044  com¬ 
puter. 

"V  =  number  of  cost  variables. 

CI  -  number  of  iterations. 


CONCLUSION 

The  Monte  Carlo  simulation  models  provide  a  tool  for  expressing  uncer¬ 
tainty  in  item  cost  estimates  and  deriving  a  corresponding  expression  of  un¬ 
certainty  in  total  costs.  The  user  must,  however,  be  aware  of  the  potential, 
the  costs,  and  the  limitations  of  the  models  before  applying  them  to  solve  a 
particular  problem. 
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WEIBULL  INPUT  REQUIREMENTS 
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The  user  of  the  Weibull  model  must  provide  four  types  of  input: 

(1 )  Cost  model  specifications 

(2)  The  number  of  cost  inputs 

(3)  An  arbitrary  6-digit  integer 

(4)  Cost  data 

Tnese  four  types  of  input  will  be  considered  in  turn. 

The  user  must  supply  the  cost  model  to  the  first  computer  program 
(WKiBl’P)  in  the  position  shown  in  Fig.  Al.  This  cost  model  is  subject  to 
several  restrictions.  First,  it  must  be  specified  in  FORTRAN  lV.7  Second,  the 
too.I  cost  must  be  called  output  (1HIST).  Third,  the  cost  variables  must  be 
called  FDATA(I)  where  1  is  any  integral  value  from  1  up  to  the  number  of  cost 
inputs  that  the  user  specifies. 

An  example  of  a  user-specified  cost  model  is  shown  in  Fig.  A2. 

The  second  and  third  types  of  input  must  be  provided  on  a  punched  card 
immediately  following  the  first  program,  WEIBTP.  The  format  of  that  card  is: 

Columns  2  to  4:  An  integer  signifying  the  number  of  cost  inputs.  This 
number  may  take  any  integral  value  from  1  to  100 

Columns  5  to  10:  An  arbitrary  6-digit  number  used  to  generate  uniform 
random  numbers 

Columns  11  to  80:  Blank 

The  final  type  of  input  provided  by  the  user  is  the  cost  data.  For  each 
cost  input  the  user  must  provide  the  most  likely  value,  a  low  value  and  the  prob¬ 
ability  that  the  actual  value  will  be  lower,  and  a  high  value  and  the  probability 
that  the  actual  value  will  be  lower.  Each  value  and  probability  is  written  on 
a  pu  ich  card  with  the  following  type  of  format: 

±X.YYE±ZZ 


where  ±X.YY  raised  to  the  power  +ZZ  is  the  value  or  probability.  For  example, 
-2700  would  be  written  as  -2.70E+03.  The  card  format  for  each  cost  input  is: 

Columns  2  to  10:  Low  value 

Columns  12  to  20:  High  value 

Columns  22  to  30:  Most  likely  value 

Columns  32  to  40:  Probability  that  actual  value  will  be  higher  than  the 
high  value 

Columns  42  to  50:  Probability  that  actual  value  will  be  lower  than  the 
low  value 

Each  cost  data  card  is  placed  in  a  sequence  such  that  its  position  in  the 
stack  of  cost  data  cards  is  equal  to  the  subscript  of  the  cost  input  variable  it 
represents.  Thus  the  card  containing  FDATA(5)  would  be  the  fifth  in  the  stack 
of  cost  data  cards  and  that  containing  FDATA(15)  would  be  fifteenth. 
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Fig.  A1 — Position  of  the  Cost  Model  in  the  Weibull  Program,  Example 
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Fig.  A2 — User-Specified  Cost  Model,  Example 
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To  further  elucidate  the  input  requirements  the  following  example  is 
given.  At  a  very  high  level  of  aggregation  the  total  recurring  costs  of  an 
infantry  division  can  be  thought  of  as  the  sum  of  operating  costs,  replacement 
costs,  and  pay. 

For  a  planned  infantry  division  the  most  likely  annual  operating  cost 
might  be  $1.0  million,  but  there  is  a  10  percent  chance  that  operating  costs 
might  be  higher  than  $2.0  million  or  lower  than  $0.8  million.  Similarly,  annual 
replacement  costs  most  likely  would  be  $0.5  million  with  a  10  percent  chance 
of  being  lower  than  $0.4  million  or  higher  than  $0.8  million.  Annual  pay  most 
likely  would  be  $0.5  million  and  would  have  a  25  percent  chance  of  being  $0.4 
million  and  a  20  percent  chance  of  being  $0.6  million. 

The  user-supplied  cost  model  is  shown  in  the  equation 

OUTPUT(IHIST)=FDATA<l)+FDATA(2)+FDATA(3) 

OUTPUT(IHIST)  is  total  recurring  cost,  FDATA(l)  is  assumed  to  be  operating 
cost,  FDATA(2)  is  replacement  cost,  and  FDATA(3)  is  military  pay. 

Next  the  user  must  supply  a  data  card  containing  the  number  of  cost 
inputs  and  a  6-digit  random  number.  An  example  of  such  a  card  is  shown  in 
Fig,  A3. 

The  3  in  col  4  signifies  that  there  are  three  cost  inputs  (operating  costs, 
replacement  costs,  and  military  pay).  The  number  987654  is  an  arbitrary 
integer  used  to  generate  random  numbers. 

Finally,  the  user  must  supply  the  input  cost  information.  Since  FDATA(l) 
refers  to  operating  costs,  the  first  input  cost  data  card  deals  with  operating 
cost.  Fig.  A4  shows  the  required  input  cost  data  card. 

Note  that  the  low,  high,  and  most  probable  values  are  written  in  the 

form: 

0.8  x  10* 

2.0  x  10s 
1.0  x  106 

The  input  cost  data  cards  for  replacement  costs  and  pay  are  similarly 
prepared.  The  user-provided  data  (excluding  the  cost  model)  for  this  example 
are  shown  in  Fig.  A5. 

To  place  the  inputs  in  proper  perspective  to  the  rest  of  the  computer 
program,  Fig.  A6  shows  the  program  deck  configuration.  The  user-provided 
inputs  are  indicated  by  arrows  and  the  7044  control  cards  are  noted  by 
asterisks. 


BETA  INPUT  REQUIREMENTS 

The  Beta  model,  like  its  Weibull  counterpart,  requires  four  types  of 

input: 

(1)  Cost  model  specifications 

(2)  The  number  of  cost  inputs 

(3)  An  arbitrary  6-digit  integer 

(4)  Cost  data 
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Fig.  A3 — Initial  Data  Card,  Weibull  Program 
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Fig.  A4 — Weibull  Input  Data,  Example 
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Fig.  A5— User-Provided  Dolo,  Weibull  Program 
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The  first  three  types  of  input  are  identical  to  those  required  for  the 
Weibull  model.  The  cost  data  requirements  are  similar  but  not  identical. 

The  user  must  imbed  within  the  first  program  deck  (called  BETATP) 
a  cost  model.  Figure  A7  shows  the  position  of  the  model  within  the  program. 

The  cost  model  is  subject  to  a  number  of  restrictions.  First,  it  must  be 
written  in  FORTRAN  IV.7  Second,  the  total  cost  must  be  called  output  (IHIST). 
Third,  the  cost  inputs  must  be  called  FDATA(l)  where  I  is  any  integral  value 
from  1  up  to  the  number  of  cost  inputs  that  the  user  specifies. 

An  example  of  a  user- specified  cost  model  is  shown  in  Fig.  A8. 

The  second  and  third  types  of  user-supplied  inputs  must  be  placed  on 
a  punched  card  immediately  following  the  third  subroutine,  EVLINT.  The 
format  of  that  card  is: 

Columns  2  to  4:  An  integer  signifying  the  number  of  cost  inputs.  This 
number  must  be  a  positive  integer  no  greater  than  100. 

Columns  5  to  10:  An  arbitrary  6-digit  number  used  to  generate 
uniform  random  numbers. 

Columns  11  to  80:  Blank. 

The  final  type  of  input  to  be  furnished  by  the  user  is  cost  data. 

For  each  cost  input  the  user  must  provide  the  most  likely  value,  a  low 
value,  and  a  high  value.  Additionally,  the  user  must  choose  a  standard  Beta 
distribution  that  best  matches  his  conception  of  how  the  probability  distribution 
would  look.  For  each  cost  input  the  above  cost  data  are  written  on  a  punched 
card  with  the  following  type  of  exponential  notation: 

±X.YYE±ZZ 

±  zz 

which  is  equivalent  to  ±X.YY  ,  For  example,  +34000  would  be  written  as 
+3.40E+04.  The  distribution  type  is  not  written  using  the  above  exponential 
notation  but  is  represented  as  an  integer. 

The  punched  card  format  for  each  cost  input  is: 

Columns  34  to  42 :  Low  value 

Columns  44  to  52 :  High  value 

Columns  54  to  62 :  Most  likely  value 

Column  63:  Distribution  type.  This  type  is  represented  by  an  integer 
with  possible  values  ranging  from  1  to  9.  Figure  A9  shows  the  type  of  dis- 
i  but  ion  corresponding  to  each  integer. 

The  cost  input  punched  cards  are  placed  immediately  after  the  punched 
card  containing  the  number  of  cost  inputs  and  the  arbitrary  integer. 

Each  cost  input  card  is  placed  in  a  sequence  corresponding  to  the  sub¬ 
script  of  the  cost  input  variable  it  represents.  For  example,  the  cost  input 
card  containing  the  costs  associated  with  FDATA(5)  would  be  fifth  in  the 
stack  of  cost  input  cards. 

To  illustrate  the  user-supplied  inputs  assume,  for  example,  that  the 
total  recurring  cost  of  an  infantry  division  is  the  sum  of  operating  cost, 
replacement  cost,  and  pay.  The  most  likely  annual  operating  cost  of  a  planned 
infantry  division  might  be  $1.0  million  and  the  upper  and  lower  limits  are  $2.0 
million  and  $0.8  million  respectively.  Since  it  is  more  probable  that  operating 
°ost,  if  not  $1.0  million,  will  be  higher  and  since  operating  cost  is  quite  volatile, 
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Fig.  A7 — Position  of  the  Cost  Model  in  the  Beta  Program,  Example 
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Fig.  A8— User-Specified  Cost  Model,  Example 
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S  FIELD 


the  cost  analyst  might  choose  distribution  type  3.  Similarly,  the  annual  re¬ 
placement  cost  most  likely,  low,  and  high  values  are  $0.5  million,  $0.4  million, 
and  $0,8  million.  Like  operating  costs,  replacement  cost  is  quite  volatile  and 
if  the  most  likely  value  is  wrong  it  is  probably  too  low.  Therefore  the  cost 
analyst  might  choose  type  3  once  again.  Finally,  the  most  likely,  low,  and  high 
values  for  military  pay  are  $0.5  million,  $0.4million,  and  $0.6  million.  Military 


Fig.  A9— Distribution  Types 


pay  does  not  change  much  and  the  most  likely  estimate  is  equally  likely  to  be 
too  high  or  too  low.  Thus  the  cost  analyst  might  consider  distribution  type  8 
to  be  the  closest  approximation  to  his  concept  of  the  military  pay  distribution. 

To  run  the  above  data  in  the  Beta  model,  the  user  must  first  supply  a 
cost  model.  The  required  cost  model  is 

OUTPUT(IHIST)=FDATA(l)+FDATA(2)+FDATA(3) 

OUTPUT(IHIST)  is  the  total  recurring  cost,  FDATA(l)  is  operating  cost, 
FDATA(2)  is  replacement  cost,  and  FDATA(3)  is  military  pay. 

Next  a  punched  card  containing  the  number  of  cost  inputs  and  a  6-digit 
random  number  must  be  supplied.  Such  a  card  is  shown  in  Fig.  A10. 

The  3  in  col  4  of  that  figure  specifies  to  the  model  that  there  are  three 
cost  inputs  (operating  cost,  replacement  cost,  and  military  pay).  The  number 
123456  is  an  arbitrary  number  used  to  generate  random  numbers. 

Finally,  the  user  must  supply  the  input  cost  Information.  FDATA(l) 
refers  to  the  first  variable  in  the  cost  model.  Examples  of  user- supplied 
input  cost  data  cards  are  shown  in  Fig.  All. 


r.st  t  -f '  w4- 


To  place  the  inputs  in  proper  perspective  to  the  rest  of  the  computer 
program,  Fig.  A12  shows  the  program  deck  configuration.  The  user-provided 
inputs  are  indicated  by  arrows  and  the  7044  control  cards  are  noted  by 
asterisks. 


ERROR  MESSAGES 

The  error  messages  for  the  Weibull  and  the  Beta  programs  are  identical 
except  for  one  message  dealing  with  input  distribution  types  for  the  Beta  pro¬ 
gram.  This  section  lists  the  error  messages  (data  errors),  explains  the  mean¬ 
ing,  and  suggests  methods  for  correcting  the  errors. 

DATA  ERROR-MODE  NOT  BETWEEN  HIGH  AND  LOW  IN  DATA  SET  XXX 

Meaning:  There  is  an  error  in  cost  data  card  XXX.  The  most  probable 
value  for  cost  variable  XXX  is  either  equal  to  or  higher  than  the  high  value 
specified  or  equal  to  or  less  than  the  low  value  specified.  All  data  are  read 
but  no  simulation  is  run  and  no  output  is  produced. 

Corrective  action:  Check  to  see  that  the  data  are  punched,  from  left  to 
right,  in  the  following  order:  low,  high,  and  most  likely.  A  second  source  of 
error  to  check  is  the  exponent  of  one  or  more  of  the  three  data  items.  Finally, 
check  to  see  that  the  data  are  punched  in  the  correct  columns. 

DATA  ERROR-LOW  GREATER  THAN  HIGH  IN  DATA  SET  XXX 

Meaning:  There  is  an  error  in  cost  data  card  XXX.  The  low  value 
specified  on  the  card  is  greater  than  or  equal  to  the  high  value.  All  data  are 
read  but  no  simulation  is  run  and  no  output  is  produced. 

Corrective  action:  Same  as  for  the  previous  error  message. 

DATA  ERROR-DISTRIBUTION  TYPE  FOR  DATA  SET  XXX 
IS  ZERO  OR  NOT  SPECIFIED 

Meaning:  This  message  is  only  generated  by  the  Beta  program.  The 
type  parameter  on  cost  data  card  XXX  is  either  not  specified  or  is  specified 
as  zero.  All  data  are  read  but  no  simulation  is  run  and  no  output  is  produced. 

Corrective  action:  Check  cost  data  card  for  the  missing  parameter. 

The  parameter  may  be  punched  in  the  wrong  column. 

ERROR  MESSAGE-RANGE  OF  TOTAL  COST  ESTIMATES  IS  LESS  THAN  TEN 
PLOT  DELETED 

Meaning:  The  simulated  total  cost  estimates  have  a  range  smaller 
than  10.0.  In  this  case  the  X  axis  cannot  be  designed  and  the  plot  is  deleted. 
Program  execution  continues. 

Corrective  action:  The  cost  data  cards  may  be  altered  by  increasing 
all  the  exponents, 
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COMPUTER  HARDWARE  AND  SOFTWARE  REQUIREMENTS 

This  section  examines  the  computer  facilities  required  to  run  the 
Monte  Carlo  programs.  Since  many  of  the  requirements  are  themselves  a 
function  of  the  kind  of  computer  used,  the  following  discussion  should  be  con¬ 
sidered  only  a  guideline. 

Hardware 

Four  hardware  features  are  required  for  the  Monte  Carlo  programs. 

First,  there  must  be  some  method  of  reading  punched  cards  since  the 
program  and  data  inputs  are  on  punched  cards. 

Second,  there  must  be  sufficient  main  (core)  storage  to  contain  the 
program.  No  precise  storage  requirements  can  be  formulated  because  of 
variations  among  computers;  however,  it  can  be  stated  that  the  variable  arrays 
employed  in  the  programs  require  3000  words  and  the  programs  themselves 
are  at  least  as  large.  Ten  thousand  words  of  storage  ought  to  be  sufficient. 

Third,  there  must  be  some  auxiliary  storage  device  to  hold  intermediate 
results.  This  device  might  be  a  tape  drive,  a  disk,  or  a  drum.  The  capacity 
of  the  storage  device  muBt  exceed  2000  words. 

Finally,  there  must  be  a  printer  to  print  the  output.  The  printer  must 
have  a  capacity  of  at  least  120  characters  per  line. 

Software 

The  programs  are  written  entirely  in  7040/7044  FORTRAN  IV.  They  do 
not  use  either  “Print"  or  “Punch*  statements. 

The  programs  are  designed  to  run  under  the  7040/  7044  operating  sys¬ 
tem.  Use  with  any  other  computer  will  entail  the  user’s  supplying  a  different 
set  of  computer  control  cards. 


PROGRAM  LOGIC 

The  purpose  of  thi  n  section  is  to  provide  sufficient  information  about 
the  programs  per  se  to  allow  them  to  be  modified  successfully.  A  verbal 
description  of  the  Weibull  program  and  subprograms  is  supplied  as  well  as  a 
description  of  how  they  are  linked  together.*  Complete  program  listings  and 
flow  charts  are  provided  for  both  programs.  The  verbal  description  is  keyed 
to  the  program  listing.  References  are  made  to  specific  lines  of  computer 
coding.  Knowledge  of  the  flow  charts,  on  the  other  hand,  is  not  specifically 
required  to  read  the  verbal  description.  However,  the  flow  charts  should 
prove  quite  useful  as  a  general  guide  to  the  programmer. 

Weibull  Program 

There  are  three  arrays  used  in  the  Weibull  program  (called  WEDBTP). 
The  first,  FDATA,  is  the  array  containing  the  values  of  the  cost  input  data. 
The  subscript  100  sets  the  maximum  number  of  cost  inputs  that  may  be  used 


♦The  Beta  program  is  similar  in  format  and  employs  identical  output  routines. 
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in  a  cost  model.  If  this  number  is  to  be  increased  then  FDATA  and  the  second 
subscript  of  WYBARY,  a  second  array  used  to  simulate  the  FDATA,  must  be 
altered  in  the  dimension  statement.  The  first  subscript  of  WYBARY  is  fixed. 
The  third  array,  output,  is  dimensioned  according  to  the  number  of  iterations 
of  the  simulations  to  be  run.  If  the  number  of  iterations  is  to  be  increased 
then  the  subscript  of  output  must  be  increased. 

There  are  six  internal  program  parameters  that  are  set  immediately 
after  the  dimension  statement.  First,  there  are  three  parameters,  101,  102, 
and  103,  that  set  the  device  numbers  for  the  system  input  device,  the  system 
output  device,  and  an  intermediate  storage  device.  The  fourth  parameter, 

ITER,  is  set  to  the  number  of  iterations  of  the  simulation  to  be  run.  If  ITER 
is  increased,  then  the  subscript  of  the  array  output  must  be  increased.  (An 
array  in  the  plot  routine  must  also  be  increased.)  The  fifth  parameter, 
NPAGES,  refers  to  the  number  of  pages  each  plot  is  to  occupy.  MI,  the  sixth 
parameter,  sets  the  number  of  class  intervals  to  be  used  in  constructing  a 
frequency  distribution  of  the  simulation  output.  If  NI  is  changed,  a  number  of 
arrays  in  the  plotting  routine  (discussed  below)  must  be  increased  also. 

After  all  internal  parameters  are  set,  two  external  parameters  are 
read:  NSETS,  the  highest  subscript  of  FDATA  actually  used  by  the  cost  model, 
and  IRANDM,  an  arbitrary  6-digit  integer  used  to  generate  random  numbers. 

After  the  above-read  statement,  the  next  block  of  coding  (through 
statement  number  122)  deals  with  reading  the  cost  input  data  and  defining  the 
Weibull  distribution  around  each  input.  First,  the  cost  data  are  read  into  the 
computer.  Then  a  check  is  made  to  see  that  the  high  cost  is  greater  than  the 
low  cost.  If  it  is  not,  an  error  message  is  written  and  a  flag  is  set  ir.  the 
parameter  ITER  by  setting  ITER  equal  to  1.  A  second  check  is  made  to  see 
if  the  mode  falls  within  the  range  of  the  high  and  low.  If  it  does  not,  an  error 
message  is  written  and  ITER  is  set  to  1.  After  the  above  checks  are  made 
ITER  is  checked  to  see  if  ITER  =  1.  If  it  does,  the  rest  of  the  data  is  read  and 
checked  but  no  simulation  takes  place.  Data  are  written  on  device  103  signaling 
that  errors  did  occur  and  the  run  is  stopped.  If  no  error  occurred,  then  initial 
parameters  for  a  Newton-Raphson  iterative  process  for  specifying  the  Weibull 
distribution  for  a  given  input  variable  are  specified.  The  iterative  process  is 
contained  in  the  coding  starting  at  statement  number  3100  and  extending  to 
statement  number  3110.  The  last  four  statements  before  statement  number  122 
transfer  the  three  parameters  specifying  the  Weibull  distribution  to  WYBARY. 
Then  control  returns  to  the  second  read  statement  and  the  process  is  repeated 
until  all  the  data  are  read  and  their  Weibull  distributions  are  specified. 

The  next  block  of  coding,  through  statement  number  9960,  is  the  simula¬ 
tion.  Basically  the  coding  consists  of  a  double  loop,  the  inner  running  from 
statement  number  9980  to  statement  number  121,  and  the  outer  from  statement 
number  9980  to  9960  The  inner  loop  generates  uniformly  distributed  random 
numbers  and  then  distributes  them  according  to  the  Weibull  distribution  asso¬ 
ciated  with  the  input  data.  The  process  is  repeated  until  all  input  data  have 
values.  In  the  outer  loop  the  cost  model  is  evaluated  on  the  basis  of  the  values 
for  the  input  data.  The  above  process  is  repeated  ITER  times,  i.e.,  until  there 
are  ITER  values  for  the  total  cost  output. 

The  next  block  of  coding,  to  statement  number  9961,  writes  the  number 
of  iterations,  the  number  of  class  intervals,  the  number  of  pages  for  each  plot, 
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the  number  of  data  inputs,  and  the  simulation  results  on  intermediate  storage 
device  103.  The  storage  device  is  backspread  so  as  to  be  in  position  to  be  read 
and  the  program  terminates. 

The  final  section  of  coding  contains  the  error  messages  and  format 
statements. 

Plot  and  Frequency  Distribution  Program 

The  plot  and  frequency  distribution  program  (called  HSTPLT)  is 
separate  from  the  Weibull  and  Beta  programs.  In  an  effort  to  reduce  core 
storage  requirements,  HSTPLT  is  loaded  on  top  of  the  preceding  Weibull  or 
Beta  programs  with  all  preliminary  results  being  saved  on  an  intermediate 
storage  device. 

The  internal  organization  of  HSTPLT  includes  a  main  program  to  plot 
the  results  of  the  simulation  and  a  subprogram  to  transform  the  simulation 
data  into  a  smooth  curve.  The  program  and  subprogram  are  considered  in 
turn. 

The  program  HSTPLT  is  conceptually  divided  into  two  parts.  The  first 
part  sets  up  the  output  histogram  and  calculates  the  mean  and  standard  devia¬ 
tion  of  the  data.  The  second  part  plots  the  output  data.  In  the  first  part  there 
are  four  arrays.  NFREQ  is  dimensioned  as  large  as  NI,  the  number  of  class 
intervals  in  the  histogram.  NFREQ  is  the  frequency  of  output  in  each  class 
interval.  XV AL  is  dimensioned  the  same  as  NFREQ  and  is  the  midpoint 
value  for  each  class  interval.  A  third  array,  SMOOTH,  contains  the  smoothed 
frequencies  corresponding  to  each  value  of  XV AL  and  is  also  dimensioned  the 
same  as  NFREQ.  The  final  array  is  output  having  a  dimension  of  ITER,  a 
parameter  that  is  set  in  the  Weibull  and  Beta  programs. 

There  are  five  internally  set  parameters.  The  first  three,  101,  102, 
and  103,  are  the  numbers  of  the  system  input  device,  the  system  output  device, 
and  an  intermediate  storage  device,  respectively.  The  fourth  parameter,  KOR, 
is  the  degree  of  the  polynomial  used  to  smooth  the  simulation  output.  The 
final  parameter,  IHOPE,  is  set  to  zero  initially.  IHOPE  counts  the  number  of 
times  HSTPLT  is  executed.  When  IHOPE  is  1,  the  program  is  operating  on 
the  unsmoothed  data.  When  IHOPE  is  2,  the  program  Is  dealing  with  smoothed 
data.  In  addition  to  the  above  parameters  there  is  a  data  statement  that  estab¬ 
lishes  four  graphic  characters  for  the  Y  axis  of  the  output  plots. 

After  the  internal  parameters  are  set,  the  program  reads  internal 
parameters  saved  on  storage  device  103.  These  parameters  are  described  in 
the  section  dealing  with  the  Weibull  and  Beta  program  logic.  If  the  parameter 
ITER  is  equal  to  1,  then  there  was  an  input  data  error  and  execution  is  halted. 
Otherwise  the  simulation  output  is  read. 

Statement  number  445  begins  the  main  loop  of  the  program.  IHOPE  is 
incremented  to  show  the  number  of  times  that  the  loop  is  being  executed.  The 
first  block  of  coding,  through  statement  number  91,  is  concerned  with  finding 
the  maximum  and  minimum  values  produced  by  the  simulation.  Immediately 
following  statement  number  91  is  a  block  of  coding  that  calculates  the  length 
for  the  plot  and  histogram  class  intervals.  Statement  90  computes  the  mid¬ 
point  of  each  class  interval.  After  statement  number  90,  through  statement 
number  92,  the  frequency  or  the  number  of  observations  falling  within  each 
class  interval  is  calculated.  The  method  is  basically  to  start  from  the  lowest 
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class  interval  and  to  test  whether  the  upper  limit  of  the  class  interval  is 
greater  than  or  equal  to  the  value  of  an  observation.  If  so,  the  number  of 
observations  within  the  class  interval  is  increased  by  1.  Frequency  calculation 
occurs  only  during  the  first  execution  of  the  main  program  loop,  since  during 
the  second  execution  smoothed  frequencies  have  been  supplied  by  the  smooth¬ 
ing  subroutine.  The  remainder  of  the  coding  for  the  first  portion  of  HSTPLT, 
through  statement  number  3200,  consists  of  a  straightforward  calculation  of 
the  arithmetic  mean  for  all  observations  and  of  the  standard  deviation. 

The  plotting  routine  is  a  conceptually  separate  part  of  HSTPLT.  How¬ 
ever,  technically  it  is  only  a  continuation  of  the  coding  for  HSTPLT.  There 
are  three  arrays.  First,  LABELX  contains  the  labels  used  for  the  X  axis  of 
the  output  plot.  A  dimension  of  16  is  given  to  LABELX,  with  LABELX(16)  set 
to  a  high  value,  although  there  are  only  15  labels  printed  owing  to  an  idio¬ 
syncrasy  in  the  coding  that  prints  the  X  axis.  LABELY  contains  the  Y  axis 
labels  and  LINE  contains  the  graphic  characters  used  to  generate  the  Y  axis. 

The  first  two  lines  of  coding  in  the  plotting  routine  determine  the 
number  of  printer  lines  per  class  interval.  The  number  of  lines  is  at  least  I. 
These  two  lines  in  conjunction  with  the  value  of  NI  control  the  number  of 
pages  the  output  plot  will  fill.  If  NI  is  greater  than  50,  then  the  number  of 
pages  is  NI/50  (rounding  upward  to  the  next  highest  integer). 

The  next  block  of  coding,  through  statement  number  450,  constructs  the 
scale  for  the  X  axis  and  fills  the  X  axis  labels  with  the  proper  values.  The 
fourth  line  of  coding  (the  line  preceding  statement  number  400)  checks  to  see 
if  the  range  of  the  values  for  the  simulation  output  is  greater  than  10.  If  it  is 
not,  control  is  transferred  to  an  error  message  routine.  The  coding  starting 
at  statement  number  400  and  ending  at  statement  number  420  constructs  the 
increments  for  the  X  axis  labels.  The  coding  from  statement  numbers  420  to 
432  establishes  the  value  of  the  first  label  and  the  remaining  coding,  through 
statement  450,  establishes  the  values  for  each  label. 

The  next  five  lines  of  coding,  through  statement  number  460,  constructs 
the  Y  axis  increments  (first  two  lines)  and  determines  the  Y  axis  label  values 
(last  three  lines). 

The  next  section,  through  two  lines  past  statement  number  483,  writes 
the  Y  axis  labels,  the  Y  axis  itself,  and  the  minimum  value  of  X. 

The  next  block  of  coding,  extending  through  statement  number  499,  prints 
the  X  axis  labels,  the  plot,  and  the  mean  value  of  X.  The  inner  loop  in  this 
coding  (DO  499  J=l,  INTWID)  is  essentially  unused  unless  the  number  of  class 
intervals  (NI)  has  been  reset  to  a  value  less  than  26.  The  plot  logic  will  net 
be  discussed  in  detail.  However,  the  coding  idiosyncrasy  with  respect  to 
LABELX  (16)  noted  above  will  be  explained.  After  the  last  label  has  been 
printed,  by  statement  number  471,  the  subscript  of  LABELX  is  increased  to 
16.  If  LABELX  (16)  is  less  than  the  minimum  value  of  X,  then  the  first  “IF" 
statement  after  statement  number  475  will  eventually  cause  statement  number 
471  to  be  executed  one  extra  time. 

After  completing  the  plot  the  maximum  value  of  X  and  the  standard 
deviation  of  X  are  written.  If  IHOPE  is  2,  i.e.,  if  this  is  the  second  time 
through  the  main  program,  execution  is  halted.  If  this  is  the  first  time,  the 
smoothing  routine  is  called.  On  return  from  the  smoothing  routine,  the 
smoothed  values  are  rounded  to  the  nearest  integer  and  control  is  transferred 
to  the  beginning  of  HSTPLT. 
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The  smoothing  subroutine  is  SHARE  subroutine  FOLYFT  and  documen¬ 
tation  is  available  through  the  SHARE  organization.  In  general,  the  subroutine 
fits  a  polynomial  to  the  unsmoothed  data  using  the  least -squares  criterion. 

The  dimension  of  the  arrays  XV AL,  NFREQ,  and  SMOOTH  must  equal 
NI,  the  number  of  class  intervals.  The  dimension  of  C,  SMYX,  and  AMEANX 
must  equal  KOR,  the  degree  of  the  polynorhia? being  used  to  smooth  the  data. 
The  dimension  of  A  must  be  KOR  squared  and  the  dimension  of  SUMX  must 
be  two  times  KOR. 

A  series  of  flow  charts,  Figs.  A  13  to  A19,  follows,  in  which  various 
routines  of  the  Monte  Carlo  simulation  are  portrayed. 


Fig.  A 13— General  Flow  Chart  of  Monto  Carlo  Routines 
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Fig.  A15— Flow  Chart  of  Sample  Routine 


Fig.  A16— Flow  Cliort  of  HISlu  Routine 


Fig.  A17 — Flow  Chart  of  Plotting  Routine 


Fig.  A 18 — Flow  Chart  of  POLYFT  Routine 
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Fig.  A19— Flow  Chart  of  WEIBTP  Routine 
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24  API  -API  +AMEANX(n  »C(I) 
CO  «  AMEANY  -APT 


SENTRY 

S1BSY5 

SCLOSE 


SRES  -  0.0 _ _ _ 

WRI TE ( 102  «29 1 _ 

FORMAT  (60H1 _ COST _ FRE 

* _ RES/43H _ INTERVAL _ 

DO  77  1-1. N1 _ 

SMOOTH ( 1 )-C0 _ 

DO  27  J-l.KOR  < _ 

SMOOTH) 1 >-SMOOTH( 1  )+C(J)»XVAL( I)»»J 
RES -FLOAT ( NFREO ( 1 ) 1 -SMOOTH ( I )  _ 

DSNORMAL IZAT ION  WITH  RESPECT  TO  XMAX  ~~ 

XVAL ( I ) -XVAL ( I  )»XMAX _ 

XVAL ( I ) -XVAL ( I l+XVALl _ 

FORMAT (  4F15.41  _ _ _ 

FREQ-FLOATI NFREO ( I  1  ) _ 

WRITE! 102. 301 XVAL ( 1 ) .FREQ. SMOOTH ( 1 ) .RES 

RETURN _ _ _ _ 

END _ 

If _ ~ _  HSTPLT _ _ 

5 _ 

■  S.SU10 


SMOOTHED 

FREQ) 


P0LY0680 : 
P0LYQ690  I 
POLY0700 
P0LY0710 
P0LY0720 
P0LY0756 
“P0T7ZT755  ' 

POLY0T50 
POLV6760" 
P0LY0770 
POLY0780 
P0LY0790 
P0LY0800 
°0LYQ810 
_f_'LY0820_ 
P0LYQ830 
POLY0840 
P0LY0850 
P0LY6B65 
POLY 08 70 
POLYoftftO 
POLY0890 
POLYCjSOO 
P0LY0910 
P0LY09S0 
POLY0930 
POLY0940 
POIY0V50 
POLY0960 
P0LYO97O 
POLY1030 


P0LY1090 


POLY 11 30 
P0LY1 14Q 


P0LY1 160 


POLY1270 

P0LY1280- 
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BZTA  PROGRAM  LISTING 


SJOB _ 009.303 t JOHNSON  »WW  t  WAIN _ X7131 _ _ 

XQPEN _ S.SU10.  REWIND _ _ _ _ _ 

SlBJOB  MCCOST  NODECK. _ _ _ _ _ _ 

8IBFTC  BET  ATP  NODECY. _ _ _ 

_ DIMENSION  OUTPUT(IOOO) _ _ _ _ _ _ 

COMMON /NONAME/NT ABLE .NSET St IRANDM.XTABLE ( 12 B ) iFDAT  AlToo ) . 

*  FMODE ( 100 ) t FRANCE  <100) tNTYPE! 100 ) 1 1  TER tNI tCUMl 9 » 128 ) 

~ _ EXTERNAL  SFX _ _ _ _ _ _ _ _ 

_ _ ITER-1000  _ _ _ 

r _ Nl-100 _ 

, _ NPAGES-2 _ _ _ _ _ 

_  I01-S _ _ _ _ _ ___ _ _ 

_ 102-6 _ _ _ _ _ ' 

_ 103-10 _ _ _ _ _  ' _ 

READtI01.102)  NSETS.IRANDH  _ _ _ _ 

102  FORMAT  (I».If) _ ‘ _ __ _ _ 

~D0  708  1-ltNSETS  _ 

'  READ  (101.103)  FLOW >FH I  OH t FMODE ( 1 ) tNTYPE  I  I ) _ 

FRANCE!  I  )«FHIGH-FLOW _ _ _ _ _ _ 

_ IF  (FRANCE!  I )  .LT.O.  )  GO  TO  9961 _ _ _ _ _ 

70 9  1FINTYPE( I). GE.l.AND tNTYPE ( I t.LE. 9J_ GO  TO  707 _ _ 

WRITE  (102.101)  1 _ 

t _ ITER-1 _ _ _ 

707  IF  (FMODE! 1 ) . LE. FHIGH.AND. FMODE ( 1 ) .GE.  FLOW)  GO  TO  708 _ 

_  WRITE  (102.104)  I _ 

_ ITER-1 _ 

_ 706  CONTINUE _ _ _ _  _ 

i _ IF  (I  TER  »EQ.  1  )  GO  TO  9960 _ 

X  GENERATES  CUMULATIVE  BETA  TABLES  FOR  NINE  BETA  EQUATIONS  ~ 

_  DO  70  _M  -  1 19~ _  ' _ ~ 

_  j-0 _ _ _ _ _ __ 

_  2 _ NTABLE  «  M _ 

_ DELTA-  0.0078125 _ 

A-  0.0 _ 

_ BB-DELTA _ _  _ _ _ _ 

_ DO  70  J  «  1.128 _ _ _ 

_ CALL  EVLINT(SFX.A.BBtFF) _ _ _ 

_ CUMIM.JI-FF _ _ 

_ XT  ABLE!  J)_-B8 _ _ _ ~ 

~  BB-BB+DELT A _ ; _ 

, _ 70 _ CONTINUE _ _ _ 

t _ LHlSIil _ _ _ _ _ 

i990a _ CALL  SAMPLE _ 

^C _ ; _ _ _ L 

JC _ SAMPLE  RETURNS  NSETS  VALUES  OF  FDATA  FOR  USE  AS  INPUT  TO  THE  COST  MODEL 

J. _ _ _ _ _ 

!  C  THE  COST  MODEL  IS  INSERTED  HERE _ _ _ 

' _ OUTPUT  (  1HIST)«FDATA(1)+FDATA(2)  +FDATA ( 3 )+FDATA ( 4) +FDATA (9)+ _ 

| _ «FDATA(6)fF0ATAI7) _ _ _ _ _ 

fc  _ _ _ _ _ _ _ _ _ _ _ _ _ 

.  c  « _ »»»»»» _ * _ *  *  »  »»*»«»»»»»»»» 

IX..  *  * _ »  ± _ *  »»»»»*»»»  »  «»»»»_»» 

lx - 

_ IF  (IHIST.GE.1TER)  GO  TO  9960 _ _ 

: _  __  I H I  ST  --.mST.t.1 . . . . . . •  _ 

_  GO  TO  9908 _  _ 

996Q  gflMILMUE _ __ 

MB11E  <iflaJL-U£H«MltflPAfiES»MSE15 _ _ _ 

_ WRITE  (10?)  ( OUTPUT (K)tX-l.ITER)  _ _  ______  _ _ _ 


I 


| 

1 

•1 


82 


_ _ 8ACKSPACE  1' _ 

BACKSPACE  10 _ _ _ _ __ _ 

_ _ _  CALL  exit _ 

996  i  WRITE  (102.100)  I _ _ _ _ _ _ 

_ _ 1TER«1 _ “ _ 

_ _ GO  TO  709  _ _ _ _ _ _ 

_ 50  FORMAT (  13. IX. 15) _ ' _ _  _ _ _ 

_ 5 1_  FORMAT (3(E10.3»1X)»I1) _ _ _ 

101  FORMAT  (46H  DATA  ERROR  -  DISTRIBUTION  TYPE  FOR  DATA  SET.I4.1X.24 

*H1S  ZERO  OR  NOT  SPECIFIED) 

103  FORMAT  (32X.3E10.2iIl) 

100  FORMAT (55H  DATA  ERROR —  LOW  GREATER  THAN  hIC-H  IN  DATA  SET  NUMBER. 
*13) 

104  FORMAT ( 56H  DAT  A  ERROR  —  MODE  NOT  BETWEEN  HIGH  AND  LOW  IN  OATA  SET 

*.  14) _ _ _ 1 _ 

END 

SIBFTC  SSSSSX  _ _ _ ~  _ 

Function  sfxixi  ~~ 

COMMON/NONAME /N TABLE .NSEtS . I RANDM.XTABLE ( 128 ) .FDATA( 100) • 

*  FmODEI 100 ( .FRANGE ( 100 ) .NTYPE! 100 ). I  TER .N I .CUMI9.128) 

C 

C  THIS  SUBROUTINE-USED  BY  PKLEQ—  DEFINES  BET“A  "EQUATION  PARAMETERS 
C 

olMf“NSj^N  con^Y&  (3.9) 

“  DATA  CONSTB  /l. 5  ,  . 5 .5. ITT.  35.  l.TSTTo . 66 »  . 57TT5T57I . I76TTT0 . 20 . 0Y2 

*  . 73. 2. ‘/s. 95 TS.  1.0. 3. 0.20. 0.4. 5.1  •  5 VT2. 5T4» 0.4»0.630..  1.5. 4. 5. 72.5 

*  r 

SF**C0N5TB  < 3 . NTaBlET* ( X**cONSTB ( 1 , NTABLE) ) *T f  lTO-X i **C OnsTB I2TNT 
»  ABLE)) _ _  _ 

RgTURN 

END _ _ _ ^ _ 

SIBFTC  SAMPLE _ _ _ _ _  ~  __  _  _ _ . _ _ _ 

_ SUBROUTINE  SAMPLE _  __  _ _ _ _ 

_C _ _ _ 

.  C  GENERATES  A  MONTE  CARLO  VALUE  FOR  EACH  INPUT  PARAMETER 
C 

.  COMMON/NONAME/NTABLE. NSETS. 1RANDMTxTABlE( 128) .FoTTAI 100). 

*  FMODE ( 1 00 ) » FRANGE  flOO ) .NTYPE ( 100 ) . I  TER .NI »CUM( 9.128 ) 

_  OA TA  JK/100043/ _ _ _ “ 

_ _ _ OATA  THQU/100000.  / _ "  _ _ ~ _ _____ 

_ REAL  MODTYP _ _ _ _ _ " _ 

_ DIMENSION  MODTYP ( 9 )  _  _ _ 

DATA  MODTYP/.75..5O..25..75».56..25»'.75.“.50..25/ 

DO  99  N*>1. NSETS  “ 

_ 1 F (NTYPE ( N ) .LT.1.0R,NTYPE(N).GT.9)G0  TO  151 _ 

■ _ I»NTYPE(N)  ' _ _ _ _ _ I _ _ _ _  _ _ 

_ S  MO  D  E  »  MO  D  T  Y  P  <  I  ) _ | _ _  _ ~ _ ~ 

_ _ M»NTYPE  (N)  _ _ _ _ _ 

_ _ GO  TO  10 _ _ _ ' _ ~~  _ 

C _ _ _ _ _ _ _ _ _ 

C  SINGLE  VALUED  INPUT _ _ _ _ 

_C _ _ _ _ _ _ _ _ _ 

151  FDATA(N)=FMODE(N) _  ~ _ _ _ 

GO  TO  99 

_ 10  L°IRANDM _ _____ _ 1 _ _ _ 

3  L«78125*L 

_ L  °L-  1 L/ JK)  *  JK _ _ _ 

IFlL-100 00 014.4.3 

4  XL«FLOAT(D/THOU 

5  fftANbM«L  _____ 

J  *  64 

00  13  K*1 .7  “  ~ 
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k'»k 


| 

I 

i 


I 


13  CONTINUE _ 

14  SAMPLX  -  XTABLEt J)  _ _ _ 

FOATAINI-FMODEINl+^RANGEIN)* (SAMPLX-SM5oEH _ 

99  CONTINUE _ _ 

100  RETURN  _ _ _ _ _  " 

END  _ _ _ _ _ . 

SlBFTC  INTEGR _ _ _ _ _ _ 

_ SUBROUTINE  EVLINTIF,  At  B,  Y  ) _ 

DIMENSION  SUB(lO).  WGT  (  5 ) 

OAT  A _ SU8/.  1 30467  360E-1,.  6746831 7QE-1. _ 

A  « 160295216 ,  .283302303,  *425562831.  .986953264. _ ‘ 

B  . 932531683.  .839704784.  .716697697,  .574437169/ 

DATA  WGT/.333356722E-1,  .747256746E-1 » 

A  .109543181.  .134633359,  .147762112/ 

DX-B-A 

Y-0.0 

00  100  I«1.5 _ 1 _ 

_ X 1-SUB (  1  )  »DX+A _ _ 

_ X2-SUB(  I+51+DX+A _ 

100  Y-Y+WGT(H«<F(Xl)+F(X2)l _ 

_ Y°Y*OX _ 

_ RETURN _ _ _ 

_ END _ 

SENTRY _ BSTATP _ _ 

_ 7664325 _ 

_ +1.00E+04  +1.60E+04  +1.20E+Q46 

+1.9QE+Q4  +2.8QE+Q4  +2 .  1QE+043 

_ _ _ _ _ +1.20E+Q4  +3.66E+Q4  +1.80E+043 

_ +1.001+04  +1.4JE+04  +1 »  10E+046 

SIBSYS _ _ _ . _ _ _ 

SIBJOB  COMPLT  NOOECX  _ _ _ _ _ 

IBFTC  HSTPLT  NODECK _ 

EXTRACT  OF  H1STO _ _ 

_ LOl-5 _ _ _ _ _ 

102-6 _ , _ _ _  _ _ . _ 

103-10 _ _ _ 

DIMENSION  NFREQI 100) ,XVAL ( 100 ) .SMOOTH! 100 > .OUTPUT < 1000  ) _ 

INTEGER  BLANK, X. DOT. APOS _ 

DATA  BLANK, X, DOT. AP0S/1H  . 1HX ♦ 1H. , 1H » / _ 

READ  < 103 )  I  TER .NI , NP AGES, N SETS _ 

IF  (ITER  . EG) .  1  )  STOP _ _ 

_ READ  ( 103 )  ( OUTPUT  (  K )  , K» 1 > I  T  ERJ _ _ _ 

_ I  TER25-I  TER/25 _ _ _ 

_ ^PR?-2.5  _ _ _ _ 

C  CALCULATES  MEAN,  STANDARD  DEVIATION.  AND  PREPARES  HISTOGRAM 

_ I HQPE-0 _ _ 

_ 445  I  HOPE- 1 HQPE  +  1 _ _ _ 

C  CALCULATE  INTERVAL  RANGE  <  R 1  NT )  FOR  N1  INTERVALS _ 

_ XM1N  -  OUTPUT! 1) _ 

_ &MASL  *  flUIEULLU _ 

_ DO  91  K  n  1,1  TER _ 

_ 1 F (OUTPUT ( X )  .LT.XMIN )  XM I N-OUTPuT ( K  ) _ 

_  1 F  (OUTPUT  (K.)  .GT.XMAX)  XMAX-OUTPuT  ( K ) _ 


<6£> 
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91  CONTINUE _ 

RINT  =  (XMAX-XMIN) /FLOAT(NI) 

CALCULATE  MEAN  VALUE  (XVAL)  FOR  EACH  INTERVAL 


XVAL  !  1 )  =XM1N  +  R1NT/2. _ _ _ 

_ DO  90  N  *  2  <N  I _ _ 

90  XVAL(N) *XVAL(N~1 >  +  RINT 

'  DETERMINE  DATA  FREQUENCY  FOR  EACH  INTERVAL 
I F ( IHOPE.EQ.2l"  60  TO  93 
DO  20  N  «l.Nl 


lltMifiarUrtljiBa 


DO  T2~  K  « 
DU“T1  n  Tl 


IM=XVAL(N)+NINI 
IF ( XLIM . GE .OUTPUT (X) )  GO  TO  92 
21  CONTINUE 


N  *  NI 


92  NFREQ(N)*NFREQ(N)+1 

C  CALCULATE  MEAN  VALUE  (XMEANI 

93  SUMX  *  0. 

DO  30  X  =  1.1  TER 
30  SUMX=SUMX+OUTPUT(X) 


aaaag BBmmaiuiaaaMMwp 


SUMSQ  *  0. 

DO  40  X  =  T7TTER 


40  SUMSQ*SUMSQ+ ( OUTPUT ( X ) -XMEAN  >** 


IF  (ITER.NE.l)  GO  TO  3210 


XSDEV-0. 


GO  TO  3200 


3210  XVAR  =  SUMSQ/ t FLOAT ( 1  TER ) - 1 • ) 


XSDEV  =  SQRT(XVAR) 


3200  CONTINUE 


DIMENSION  LABELXU5)  .LA8ELY  <  10 )  .L INE  ( 132  ) 


DETERMINES  WIDTH  OF  PLOTTING  INTERVAL 


I NTW I D=50/N  I _ 

I F (  INTWID.EQ.O)  INTWID=1 


XRANGEnXMAX-XMIN 


I FACTR=0 


h  nwmmtiaiiw 


IFIXINCR.LT. 1. )  GO  TO  498 


400  I F ( XINCR.GE.1..AND.XINCR.LE.10. )  GO  TO  410 


XI  NCR* XI  NCR/10. 


IFACTR=IFACTR+1 


410  I NCRX*X I  NCR 


IF(  INCRX.E0.3)  INCRX-2 
I F  <  INCRX.GTT3 . ANDTTNCRX.LE.7 )  INCRX*5 


1 F (  INCRX.EQ. 

T~5C 


IIIIN  H  1 1 1 1 1 WM 

1MDI1UWI: 


1 NCRX* 10*1 NCRX 


420  CONTINU 


SCALOW-SCALOW/IO. 


IFACTR=IFACTR+1 

GO  TO  425 _ 

430  LOSCAL-SCALOW 

IF(  IFACTR.EQ.O)  GO  TO  440“ 

DO  440  I =  1 » I FACTR 
LO5CAL*10*LOSCAL 


440  CONTINUE 

431  I F (FLOAT (LOSCAL+ 1NCRX )  ,GE«  XMIN)  GO  TO  432 _ 

LOSCAL-LOSCAL+I NCRX 

GO  TO  431 

432  CONTINUE  _ _  J 

LABELX U)«LOSCAL+INCRX _ '  _ _ 

DO  450  1=2*15 


450  LABELX  (  I  )=LABELX<  I-1J+INCRX 


LABELX  NOW  CONTAINS  15  INCREMENTAL  VALUES  TO  BE  USED 
FOR'  SCALE  MARKINGS  ALONG  THE  X  AXIS 


DESIGN  Y  AXIS  SCAL E 

M AX YM  I  T E R » S  )  / Nl 

I NCRY=MAXV/ 10  ' _ 

LABELYI 1)=INCRY 


FOR  SCALE  MARKINGS  ALONG  THE  Y  AXIS 


PRINT  Y  SCALE 


WRITE  ( 102 . 402 ) < LABELY <  I  > . I -1 . 10 ) 


DO  469  1-16*116.10 


469  L I NE  <  I  )  “DOT 


DO  481  1-1.15 


481  LINEUI-BLANK _ 


DO  482  1-17.115*2 


482  L I Nf ( I ) -BLANK _ 

_  DO  483  1ST  ART-18. 108.10 _ _ 

_ 1  STOP- 1  STAR T+6 _ ' 

DO  483  I -1ST ART . I  STOP  .2 


483-UM.I.IJ.lAgQ 


itimirnwpirnri 


BEGIN  PLOTT I NG 


J.ABELING  X  AXIS 


DX  -  XRANGE/  ( FLOAT (I NTWI D ) «FLOAT ( N I ) ) 
DY  -  FLOAT(MAXY)/100 


XVALU-XMI 


AB-1 


■  Elk'll 


DO  499  J-1.INTWID 


YVALU-0. _ _ _ 

_  1-1 _ _ 

477  IF  <FLOAT<NFREQ<KJ I.LE.YVALU)  GO  TO  475' 


L1NE(I)-BLANK _ _ _ _ 

i  “i*i _ _ _ _ _ _  i 

YVALU=YVAlU*-DY 
GO  TO  477 

475  LINE ( I ) =X _ 

XVAL.U=XVAlU+DX 


c  print  x  value  If  Applicable 
c 

- 1 F  (FLOAT-! LABEL X  ( LAB)  )  .GTV  XVAIILH 


c  PRINT  xmean  IF  IT  OCCURS  in  this  interval 


IFIXMEAN.lt. XVALU  .OR.  XMEAN. GE .< XVALU+DX ))  GO  TO  471 

GO  TO  472 _ _ _ _ _ 

470  IF  (XMEAN.LT. XVALU  .ORw  XMEAN. 6fc.  ( XVALU+DX ) )  GO  TO  473 

GO  TO  474 _ _ _ 

471  WRITE! 102. 405)LABELX( L AB ) , ( L INE  1 1 1 )  ,  1 1  -1 , 1  ) _ 

LAB-LAB+1 _ __ _ _ _ .  _ 

GO  TO  499 

472  WRITE! 102 ,406 ) LABELX ( LAB ) .XMEAN .< LI NE< II). 1 1>12. 1) 

LAB-LAB+1 

GO  TO  499 

473  WRITE  (102.407)  ( L INE ( 1 1 ) . 1 1 = 1 . 1 ) _ _ 

_ GO  TO  499 _ ' _ " _ _ _ _ 

474  WR1TE(I02.408)XMEAN.(LINE(II),II=12,I) _ . 

499  CONTINUE _ _ _ ~  _ 

_ WRITE!  T 02.411) XMAX _ _  _ _ _ _ 

_ WR ITEI  I02 ,41 2  )  XSPEV _ _ ' _ _ _ _ 

412  FORMAT  (Z//.9H  STD  DEV-.F10.3) _ _ _ 

GO  TO  99 _ _ _ 

I  498  WRITE!  102.409) _ _ __ _ '  _ _ _ 

99  CONTINUE _ 

I  F  (  IH(5pe»NE.  1 )  CALL  EXIT _ _  _ _ _ 

"  __  CALL  POLYFT  (XVAL. NFREQ, NI .KOR. SMOOTH) _ _ _ 

DO  505  I  *1  ,N  I _ _ _ _ ' _ _ _ 

505  NFREQ (  I  )  “SMOOTH  (  I )  +  .500 _ ' _ _ 

_ GO  TO  445 _ _ _ _  “ _  _ _ _ 

402  FORMAT (1H1.12X.4H  OOP . 10 ( 6X  .  14  )  )  _ 

403  FORMAT  (116A1/) _ _  _ _ _ 

_  404  FORMAT ( 6HXMIN  -.F10.3//) _ _ _ _ _ 

405  FORMAT!  16. 9X.100A1  ) _ _ _ _  " _ 

_  406  FORMAT ( I 6 ,3X .7HXMEAN  - . F 1 0 .3.88A1 )  _ 

407  FORMAT! 15X.100A1 ) 

408  FORMAT (9X.7HXMEAN  ■ .F1Q.3.88A1 )  _ 

409  FORMAT  ( 79HOERROR  MESSAGE  -  RANGE  OF  TOTAL  COST  ESTIMATES  IS  LES 

_ »S  THAN  TEN.  PLOi  DELETED  )  '  _ | _ 

411  FORMAT ( 6HXMAX  -^10.3) _ _ ' _ _ 

END. _ _ _ _ _ _ 

SIBFTC  FTKURV _ _ _ _ _ 

_ SUBROUTINE  POLYFT ( XVAL .NFREQ ,N I .XOR, SMOOTH) _ 

. _ I01«5 _ _ _ _ 

_ 102-6 _ _ _ _ _ 

_ 103-10 _ ' _ 


C  POLY 00 10 

C  DIMENSION  FOR  ARGUMENTS  ~  POLY0020 

DIMENSION  XVAL(100).NFREQ(100),7i!625).C< 25). SMOOTH! 100) _ 

C _ _ P0LYQQ40 

C  DIMENSION  FOR  SELF-GENERlTfED  VALUES  POLY0050 

DIMENSION  SUMimoO)  .SMYX!  50).  AMEANX!  50)  POLYOO69 
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NORMALIZATION  WITH  RESPECT  TO  XMAX 

XVALl-XVAL(l) _ _ 

DO  707  NOJ-l.NI _ 

XVAL ( NO J ) -XVAL ( NOJ I -XVAL 1 _ 

XMAX-XVALINI ) 


DO  102  I *1 »N I 


XVAL  ( 1 1 “XVAL ( I) /XMAX 


P0LY0190 


F ORMULATION  OF  NORMAL  EQUATIONS _ 

DO  3  J^_l  *X TOR 

SUMX ( J ) -SUMX ( J ) +XVAL ( 1 1 »*J _ 

DO  A  I-1»NI _ _ _ _ 

SUMY-SUMY+FLOAT (NFREQ( I ) ) 


AMEANY-  SUMY/PTS 


DO  6  J-l. XOR 


AMEANXI J ) ■  SUMX  (J)/PT$ _ 

D06  1-1 *NI 
SMYX  ( J ) -SMYX  <  J I +FL0AT ( NFREQ (II) *XVAL< 1 1 **J 


DO  __8_I « 1  »XOR _ _ 

Cl  II-  SMYX ( I )  -PTS»AMEANX 1 1 ) »AMEANY 
DO  8  J-l, XOR  _ 


X-  I+J 


IJ  -  C  J-l  )»XOR  +  1 _ 

A ( I J i  -SUMX ( X )  -PTS*  AME ANX ( I ) » AMEANX 1 J ) 


POLY0270 

POLY0280 

POLY0290 

POLY0300 


CROUT ♦  i 


DEDUCTION  METHOD 


1 1 1 - ( 1-1 )»XQR  +  1 _ 

11  A 1 1 1 1 )  -  AC  1 1 1  1 /  A<1) _ 

DO  12  J-2  >X0R _ 

XM-  J-l 


DO  14  I-J.XOR 


_ API-  0.0 _ 

_ DPI 14  X-l.XM _ _ 

IK  -IK-1 )«  XOR+  J _ 

_ XJ  «(J-1)»  XOR+  X _ 

114  API  -  API  +  A ( 1 X I  *A(KJ) 


_ I  J  - (  J-l )«  XOR+  1 _ 

(IJ)  -  A  ( 1 J I  -APJ 


P-  J+ 


44*  443 


_  XI  -d-l)»  XOR  +  X _ 

116  API  -API  +A(JX)  *A (XI ) 


Jl  »I1-1H  XOR  ±  J _ 


POLY0350 


P0LY0360 


POLY0370 


PQLY0400 
PQLY0410 
PQLY0420 
PQLY04I5 

_ P0LYQ44Q 

POLYQ45d 

_ PQLY0460 

PQLYQ470 

_  PQLY0480 

~P0LYQ49Q 
PQLY0500 
P0LY0510 
POLY0520 


P0LY0530 


P0LY0540 
P0LYQ550 
PQLY0360 
P0LYQ570 
POLY0580 


_ POLY0590. 

PQLY060 


POLY0610 


P0LY062 


POLYO630 


P0LY0640 

PQLY063Q: 


PQLY0670: 

POLY0680 


PQLY0690 


JJ  .(J-mnc5R  +  J 


16  A ( J 1  )  =  ( A ( JI )  -API ) /A ( J J) 


445  D  '1MY*  0.0 


P0LY0700 


P0LY0710 


POLYO720 

POLY0730 


12  CONTINUE 


P0LY0?40 


C<1)  *CI1)/A !  1 ) 
DO  18  I«2»KOR  " 


API*  0.0 


I M* I -1 


POLY0750 


POLY0760 


P0LY0770 


DO  118  K*1,IM 

POLY0780 

IK  * (K-l )*  +  t 

P6LY0790 

118  API* API  +A ! I K)  *C(K) 

POLY0800 

II  •( 1-1 )*  KOR  +  I 

POLY0810 

is  cm  -icm-  APD  /  xrnr- 

POLYQ«j20 

KORM*  KOR-1 

POLY0830 

IF  (KORM)  122*  123.122 

POLY0840 

122  DO  211*1 »k6Rm 

f>OLY0850 

API*  0.0 

POLY0860 

M*  KOR-I 

POLY0870 

MP*  M+l 

POLY0860 

DO  121  k*mp.k£5r 

POLY0890 

MK  «TK-11*  KOR  *-M 

POLY0900 

121  API  “API  +  AtMkl*  C ( K ) 

POLY0910 

21  CTmF~  =CTMT"^PT 

POLY0920 

123  API*  0.0 

POLY0930 

DO  24  1=1. KOR 

POLY0940 

24  API  =AP1  +AMEANX ( I )  *C!I> 

PCLY0950 

CO  -  AMEANY  -API 

POLY0960 

C 

POLY0970. 

778  SRES  *  0.0 

POLY1030 

WRITE! 102,29) 

29  FORMAT  (60H1  COST  FREQ 

SMOOTHED 

*  RES/43H  INTERVAL 

FREQ) 

DO  77  1*1, NI 

SMOOTH (I )*C0 

DO  27  J*1 .KOR 

POLY  1090 

27  SMOOTH! I > .SMOOTH! I ) +C ( J) *XVAL ( I ) **J 

RES-FLOAT 1 NFREQ ( I) ) -SMOOTH ( I ) 

C 

POLY 11 30 

c  DENORMALIZATION  with  respect  to  xmax 

POLY 11 40 

XVAL  < I ) =XVAL ( I ) *XMAX 

XVAL 1 I ) *XVAL 1 I l+XVALl 

30  FORMAT!  4F15.4) 

POLY  1160 

FREQ*FLOAT ( NFREQ! I )  ) 

77  WRITE!  102,30) XVAL ( I ) .FREQ, SMOOTH! I ) .RES 

RETURN 

POLY1270 

END 

POLY1280 

SENTRY  HSTPLT 

SCL05E _  S.SU10 _ 
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